home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / AA_51.ZIP / AA.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  9KB  |  473 lines
  1. /* This program calculates orbits of planetary bodies and reduces
  2.  * the coordinates of planets or stars to geocentric and topocentric
  3.  * place.  An effort has been made to use rigorous methods throughout.
  4.  *
  5.  * References to AA page numbers are to The Astronomical Almanac, 1986
  6.  * published by the U.S. Government Printing Office.
  7.  *
  8.  * The program uses as a default the orbital perturbations given in
  9.  * Jean Meeus, "Astronomical Formulae for Calculators", 3rd ed.,
  10.  * Willmann-Bell, Inc., 1985.
  11.  *
  12.  * Warning! Your atan2() function may not work the same as the one
  13.  * assumed by this program.
  14.  * atan2(x,y) computes atan(y/x), result between 0 and 2pi.
  15.  *
  16.  * S. L. Moshier, November, 1987
  17.  */
  18.  
  19.  
  20.  
  21. /* Conversion factors between degrees and radians */
  22. double DTR = 1.7453292519943295769e-2;
  23. double RTD = 5.7295779513082320877e1;
  24. double RTS = 2.0626480624709635516e5; /* arc seconds per radian */
  25. double STR = 4.8481368110953599359e-6; /* radians per arc second */
  26. double PI = 3.14159265358979323846;
  27.  
  28. /* Standard epochs.  Note Julian epochs (J) are measured in
  29.  * years of 365.25 days.
  30.  */
  31. double J2000 = 2451545.0;    /* 2000 January 1.5 */
  32. double B1950 = 2433282.423;    /* 1950 January 0.923 Besselian epoch */
  33. double J1900 = 2415020.0;    /* 1900 January 0, 12h UT */
  34.  
  35. /* approximate motion of right ascension and declination
  36.  * of object, in radians per day
  37.  */
  38. double dradt = 0.0;
  39. double ddecdt = 0.0;
  40. /* Data structures containing orbital elements of
  41.  * objects that orbit the sun.  See kep.h for the definition.
  42.  */
  43. #include "kep.h"
  44.  
  45. /* Space for star description read from a disc file.
  46.  */
  47. struct star fstar;
  48.  
  49. /* Space for orbit read from a disc file. Entering 99 for the
  50.  * planet number yields a prompt for a file name containg ASCII
  51.  * strings specifying the elements.
  52.  */
  53. struct orbit forbit;
  54.  
  55. /* Orbits for each planet.  The indicated orbital elements are
  56.  * not actually used, since they are now calculated
  57.  * from the Meeus expansions.  Magnitude and semidiameter
  58.  * are still used.
  59.  */
  60.  /* Programs to compute perturbations. */
  61. #if !BandS
  62. int omercury();
  63. #endif
  64. int cmercury();
  65. struct orbit mercury = {
  66. "Mercury        ",
  67. 2446800.5, /* January 5.0, 1987 */
  68. 7.0048,
  69. 48.177,
  70. 29.074,
  71. 0.387098,
  72. 4.09236,
  73. 0.205628,
  74. 198.7199,
  75. 2446800.5,
  76. -0.42,
  77. 3.36,
  78. #if BandS
  79. 0,
  80. #else
  81. omercury,
  82. #endif
  83. cmercury,
  84. 0.0,
  85. 0.0,
  86. 0.0
  87. };
  88.  
  89. #if !BandS
  90. int ovenus();
  91. #endif
  92. int cvenus();
  93. struct orbit venus = {
  94. "Venus          ",
  95. 2446800.5,
  96. 3.3946,
  97. 76.561,
  98. 54.889,
  99. 0.723329,
  100. 1.60214,
  101. 0.006757,
  102. 9.0369,
  103. 2446800.5,
  104. /* Note the calculated apparent visual magnitude for Venus
  105.  * is not very accurate.
  106.  */
  107. -4.40,
  108. 8.34,
  109. #if BandS
  110. 0,
  111. #else
  112. ovenus,
  113. #endif
  114. cvenus,
  115. 0.0,
  116. 0.0,
  117. 0.0
  118. };
  119.  
  120. /* Meeus purturbations are invoked for earth if BandS = 0.
  121.  * If BandS = 1, the more accurate B&S expansion is computed.
  122.  * Fixed numerical values will be used if read in from a file
  123.  * named earth.orb.
  124.  * See kfiles.c, kep.h, oearth.c, and pearth.c.
  125.  */
  126. #if !BandS
  127. int oearth();
  128. #endif
  129. int cearth();
  130. struct orbit earth = {
  131. "Earth          ",
  132. 2446800.5,
  133. 0.0,
  134. 0.0,
  135. 102.884,
  136. 0.999999,
  137. 0.985611,
  138. 0.016713,
  139. 1.1791,
  140. 2446800.5,
  141. -3.86,
  142. 0.0,
  143. #if BandS
  144. 0,
  145. #else
  146. oearth,
  147. #endif
  148. cearth,
  149. 0.0,
  150. 0.0,
  151. 0.0
  152. };
  153. extern struct orbit earth;
  154.  
  155.  
  156. int omars();
  157. int cmars();
  158. struct orbit mars = {
  159. "Mars           ",
  160. 2446800.5,
  161. 1.8498,
  162. 49.457,
  163. 286.343,
  164. 1.523710,
  165. 0.524023,
  166. 0.093472,
  167. 53.1893,
  168. 2446800.5,
  169. -1.52,
  170. 4.68,
  171. omars,
  172. cmars,
  173. 0.0,
  174. 0.0,
  175. 0.0
  176. };
  177.  
  178. #if !BandS
  179. int ojupiter();
  180. #endif
  181. int cjupiter();
  182. struct orbit jupiter = {
  183. "Jupiter        ",
  184. 2446800.5,
  185. 1.3051,
  186. 100.358,
  187. 275.129,
  188. 5.20265,
  189. 0.0830948,
  190. 0.048100,
  191. 344.5086,
  192. 2446800.5,
  193. -9.40,
  194. 98.44,
  195. #if BandS
  196. 0,
  197. cjupiter,
  198. #else
  199. ojupiter,
  200. 0,
  201. #endif
  202. 0.0,
  203. 0.0,
  204. 0.0
  205. };
  206.  
  207. int osaturn(), csaturn();
  208. struct orbit saturn = {
  209. "Saturn         ",
  210. 2446800.5,
  211. 2.4858,
  212. 113.555,
  213. 337.969,
  214. 9.54050,
  215. 0.0334510,
  216. 0.052786,
  217. 159.6327,
  218. 2446800.5,
  219. -8.88,
  220. 82.73,
  221. osaturn,
  222. csaturn,
  223. 0.0,
  224. 0.0,
  225. 0.0
  226. };
  227.  
  228. int ouranus(), curanus();
  229. struct orbit uranus = {
  230. "Uranus         ",
  231. 2446800.5,
  232. 0.7738,
  233. 73.994,
  234. 98.746,
  235. 19.2233,
  236. 0.0116943,
  237. 0.045682,
  238. 84.8516,
  239. 2446800.5,
  240. -7.19,
  241. 35.02,
  242. ouranus,
  243. curanus,
  244. 0.0,
  245. 0.0,
  246. 0.0
  247. };
  248.  
  249. int oneptune(), cneptune();
  250. struct orbit neptune = {
  251. "Neptune        ",
  252. 2446800.5,
  253. 1.7697,
  254. 131.677,
  255. 250.623,
  256. 30.1631,
  257. 0.00594978,
  258. 0.009019,
  259. 254.2568,
  260. 2446800.5,
  261. -6.87,
  262. 33.50,
  263. oneptune,
  264. cneptune,
  265. 0.0,
  266. 0.0,
  267. 0.0
  268. };
  269.  
  270. /* Note there are no perturbation formulas for Pluto.
  271.  * The program automatically uses the given numerical
  272.  * values since the pointers to perturbation subroutines
  273.  * are null.
  274.  * For some reason the J2000 orbit given for Pluto in the AA does
  275.  * not give the same results as the "of date" orbit.  Yet, results
  276.  * are the same for the other planets.
  277.  */
  278. struct orbit pluto = {
  279. "Pluto          ",
  280. 2446640.5,
  281. 17.1346,
  282. 110.204,
  283. 114.21,
  284. 39.4633,
  285. 0.00397570,
  286. 0.248662,
  287. 355.0554,
  288. 2446640.5,
  289. -1.0,
  290. 2.07,
  291. 0,
  292. 0,
  293. 0.0,
  294. 0.0,
  295. 0.0
  296. };
  297.  
  298. /*
  299. int otest(), ctest();
  300. */
  301. struct orbit test = {
  302. "Test orbit     ",
  303. 2446800.5,
  304. 1.8498,
  305. 49.457,
  306. 286.343,
  307. 1.523710,
  308. 0.524023,
  309. 0.093472,
  310. 53.1893,
  311. 2446800.5,
  312. -1.52,
  313. 4.68,
  314. /*
  315. otest,
  316. ctest,
  317. */
  318. 0,0,
  319. 0.0,
  320. 0.0,
  321. 0.0
  322. };
  323.  
  324.  
  325. /* coordinates of object
  326.  */
  327. int objnum = 0;    /* I.D. number of object */
  328. static double robject[3] = {0.0}; /* position */
  329. /* ecliptic polar coordinates:
  330.  * longitude, latitude in radians
  331.  * radius in au
  332.  */
  333. double obpolar[3];
  334.  
  335. /* coordinates of Earth
  336.  */
  337. /* Heliocentric rectangular equatorial position
  338.  * of the earth at time TDT re equinox J2000
  339.  */
  340. double rearth[3];
  341. /* Corresponding polar coordinates of earth:
  342.  * longitude and latitude in radians, radius in au
  343.  */
  344. double eapolar[3];
  345.  
  346. /* Julian date of ephemeris
  347.  */
  348. double JD;
  349. double TDT;
  350. double UT;
  351.  
  352. /* flag = 0 if TDT assumed = UT,
  353.  *      = 1 if input time is TDT,
  354.  *      = 2 if input time is UT.
  355.  */
  356. int jdflag = 0;
  357.  
  358. /* correction vector, saved for display  */
  359. double dp[3];
  360.  
  361. /* display formats for printf()
  362.  */
  363. extern char *intfmt, *dblfmt;
  364.  
  365. /* display enable flag
  366.  */
  367. int prtflg = 1;
  368.  
  369. /* Tabulation parameters
  370.  */
  371. static double djd = 1.0;
  372. static int ntab = 1;
  373.  
  374. /* Main program starts here.
  375.  */
  376. void main()
  377. {
  378. int i;
  379.  
  380. struct orbit *elobject;    /* pointer to orbital elements of object */
  381. double getdate(), gethms();
  382.  
  383. kinit();
  384.  
  385. loop:
  386.  
  387. prtflg = 1;
  388. printf( "Enter starting date of tabulation\n" );
  389. JD = getdate(); /* date */
  390. JD += gethms();    /* time of day */
  391. update(); /* find UT and ET */
  392. printf( "Julian day %.7lf\n", JD );
  393.  
  394. getnum( "Enter interval between tabulations in days", &djd, dblfmt );
  395. getnum( "Number of tabulations to display", &ntab, intfmt );
  396. if( ntab <= 0 )
  397.     ntab = 1;
  398.  
  399. loop1:
  400. getnum( "Planet number 0-9 or 88 to read star, 99 to read orbit",
  401.     &objnum, intfmt );
  402.  
  403. switch(objnum)
  404.     {
  405.     case -1: exit(0);
  406.     case 0: elobject = 0;
  407.         printf( "\n                   The Sun\n" );
  408.         break;
  409.     case 1: elobject = &mercury; break;
  410.     case 2: elobject = &venus; break;
  411.     case 3: elobject = 0;
  412.         printf( "\n                   The Moon\n" );
  413.         break;
  414.     case 4: elobject = &mars; break;
  415.     case 5: elobject = &jupiter; break;
  416.     case 6: elobject = &saturn; break;
  417.     case 7: elobject = &uranus; break;
  418.     case 8: elobject = &neptune; break;
  419.     case 9: elobject = &pluto; break;
  420.     case 10: elobject = &test; break;
  421.     case 88:
  422. morstar:    elobject = (struct orbit *)&fstar;
  423.         i = getstar( elobject );
  424.         if( i == 1 )
  425.             goto loop1;
  426.         if( i == 0 )
  427.             break;
  428.         goto operr;
  429.     case 99:
  430.         elobject = &forbit;
  431.         i = getorbit( elobject );
  432.         if( i == 1 )
  433.             goto loop1;
  434.         if( i == 0 )
  435.             break;
  436.     default:
  437. operr:        printf( "Operator error.\n" );
  438.         goto loop;
  439.     }
  440.  
  441. if( elobject == (struct orbit *)&fstar )
  442.     showcname( &elobject->obname[0] );
  443. else if( elobject )
  444.     printf( "\n                  %s\n", &elobject->obname[0] );
  445.  
  446. for( i=0; i<ntab; i++ )
  447.     {
  448. /* print Julian date */
  449.     printf( "\nJD %.2f,  ", JD );
  450.     update();
  451.  
  452. /* Always calculate heliocentric position of the earth */
  453.     kepler( TDT, &earth, rearth, eapolar );
  454.  
  455.     switch( objnum )
  456.         {
  457.         case 0: dosun(); break;
  458.         case 3: domoon(TDT); break;
  459.         case 88: rstar( elobject ); goto morstar;
  460.         default:
  461. /* calculate heliocentric position of the object */
  462.             kepler( TDT, elobject, robject, obpolar );
  463. /* apply correction factors and print apparent place */
  464.             reduce( elobject, robject, rearth );
  465.             break;
  466.         }
  467.     printf( "\n" );
  468.     JD += djd;
  469.     }
  470. goto loop;
  471. }
  472.  
  473.